#ifndef GEOMETRY_H_
#define GEOMETRY_H_

#include <iostream>
#include <vector>
#include <boost/numeric/interval.hpp>
#include <CGAL/Gmpq.h>
typedef CGAL::Gmpq gmp;
typedef boost::numeric::interval<double> interval;

class Point
{
public:
    double x, y;
    Point(double x_, double y_): x(x_), y(y_) { };
    Point() {};
};

std::istream& operator>>(std::istream& in, Point& a)
{
   in >> a.x >> a.y;
   return in;
}

class Segm
{
public:
    Point a, b;
    Segm(Point a_, Point b_): a(a_), b(b_) {};
    Segm()
    {
       this->a.x = 0;      this->b.x = 0;
       this->a.y = 0;      this->b.y = 0;
    };
};

int get_predicate(double& a, double& s, interval& b, gmp& c)
{
   if (-a > s)
        return 1;
    else 
       if(a > s)
          return 1;
    if(b.lower() > 0)
        return 1;
    else 
       if(b.upper() < 0)
         return -1;   
    if (c > 0)
        return 1;
    else 
       if (c < 0)
         return -1;

    return 0;
}

int left_turn(Point& a, Point& b, Point& c)
{
    double eps = std::abs((b.x - a.x) * (c.y - a.y)) + std::abs((b.y - a.y) * (c.x - a.x));
    double s = 4 * eps * std::numeric_limits<double>::epsilon();
    interval a_xi(a.x),  b_xi(b.x), c_xi(c.x), a_yi(a.y), b_yi(b.y), c_yi(c.y);
    gmp a_xg(a.x),  b_xg(b.x), c_xg(c.x), a_yg(a.y), b_yg(b.y), c_yg(c.y);

    double res_simple = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
    interval res_int = (b_xi - a_xi) * (c_yi - a_yi) - (b_yi - a_yi) * (c_xi - a_xi);
    gmp res_gmp = (b_xg - a_xg) * (c_yg - a_yg) - (b_yg - a_yg) * (c_xg - a_xg);

    return get_predicate(res_simple, eps, res_int, res_gmp);
};

template<class T> T pow_2(T x) 
{
   return x * x;
};

template<class T> T distance_to_segm_2(Segm& s, Point& a) 
{
   T s1_x(s.a.x), s1_y(s.a.y), s2_x(s.b.x), s2_y(s.b.y), a_x(a.x), a_y(a.y);
   return (s1_y - a_y) * (s1_x - s2_x) - (s1_x - a_x) * (s1_y - s2_y);
};

int closer_to_segm(Segm& ss, Point& a, Point& b)
{
    double eps = pow_2<double>(std::abs((ss.a.y - a.y) * (ss.a.x - ss.b.x)) + std::abs((ss.a.x - a.x) * (ss.a.y - ss.b.y)))
                          + pow_2<double>(std::abs((ss.a.y - b.y) * (ss.a.x - ss.b.x)) + std::abs((ss.a.x - b.x) * (ss.a.y - ss.b.y)));
    double s = 16 * eps * std::numeric_limits<double>::epsilon();

    double res_simple = distance_to_segm_2<double>(ss, a) - distance_to_segm_2<double>(ss, b);
    interval res_int = distance_to_segm_2<interval>(ss, a) - distance_to_segm_2<interval>(ss, b);
    gmp res_gmp = distance_to_segm_2<gmp>(ss, a) - distance_to_segm_2<gmp>(ss, b);

    return get_predicate(res_simple, eps, res_int, res_gmp);
};

template<class T> T lenght_2(Segm& s)
{
   T s1_x(s.a.x), s1_y(s.a.y), s2_x(s.b.x), s2_y(s.b.y);
   return pow_2(s1_x - s2_x) + pow_2(s1_y - s2_y);
};

int longer(Segm& s1, Segm& s2)
{
    double res_simple = lenght_2<double>(s1) - lenght_2<double>(s2);
    double eps = pow_2(s1.a.x - s1.b.x) + pow_2(s1.a.y - s1.b.y)
                          + pow_2(s2.a.x - s2.b.x) + pow_2(s2.a.y - s2.b.y);
    double filter = 16 * eps * std::numeric_limits<double>::epsilon();

    interval res_int = lenght_2<interval>(s1) - lenght_2<interval>(s2);
    gmp res_gmp = lenght_2<gmp>(s1) - lenght_2<gmp>(s2);

    return get_predicate(res_simple, eps, res_int, res_gmp);
};

int the_furthest_index(std::vector<Point>& point_list, int& i, int& fur)
{
   int size = point_list.size();
   Segm s(point_list[i], point_list[(i + 1) % size]);
   while (closer_to_segm(s, point_list[fur], point_list[(fur + 1) % size]) == -1)
   {fur = (fur + 1) % size;}

   return fur;
};

Segm find_the_furthest(std::vector<Point>& point_list)
{
    int size = point_list.size(), fur = 1;
    Segm res;
    for (int i = 0; i < size; i++)
    {
        fur = the_furthest_index(point_list, i, fur);
        Segm s1(point_list[i], point_list[fur]);
        Segm s2(point_list[(i+1) % size], point_list[fur]);
        if (longer(s1, res) == 1)
            res = s1;
        if (longer(s2, res) == 1)
            res = s2;
    }
    return res;
}
#endif